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SUMMARY 


The development of a numerical simulation of time-dependent, turbulent, 
compressible flow about two-dimensional multi-element airfoils of arbitrary 
shape is described. The basis of this simulation is a technique of automatic 
numerical generation of coordinate systems fitted to the multiple bodies re- 
gardless of their number or shape. 


Incompressible solutions have been run for NACA airfoils at Reynolds 
numbers as high as 10 6 , but present drag predictions are about twice the 
experimental values. Procedures have been developed whereby the coordinate 
lines are automatically concentrated in the boundary layer at any Reynolds 
number. The compressible turbulent solution involves an algebraic eddy 
viscosity turbulence model. The laminar version has been run for transonic 
flow at free-stream Mach numbers up to 0.9. 


INTRODUCTION 


The overall purpose of the present research project is to develop a 
numerical simulation of time-dependent, turbulent, compressible flow about 
two-dimensional multi-element airfoils of arbitrary shape. The basis of this 
simulation is a technique of automatic numerical generation of coordinate 
systems fitted to the multiple bodies regardless of their number or shape. 

This procedure eliminates the shape of the bodies as a complicating factor 
and allows the flow about arbitrary bodies to be treated essentially as 
easily as that about simple bodies. All computation can be done on a rec- 
tangular transformed field with a square mesh regardless of the shape or 
number of bodies in the physical field. 

In the effort thus far, numerical solutions have been developed for the 
time-dependent incompressible Navier-Stokes equations for laminar flow about 
arbitrary multiple airfoils, and the compressible Navier-Stokes equations 
for laminar and turbulent flow about single airfoils. Continuous refinement 
of the techniques is being made as higher Reynolds numbers are considered. 

A number of papers have reported the work so far in references 1-8, and the 
coordinate system code has been made available for general use as described 
in reference 8. 
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The development of the compressible turbulent solution is now being ex- 

* lrf °fJ 8 U8ing a new coordinate system configuration having 
the infinity boundary collapsed to a point. Refinement of this solution and 
the incompressible solution for multiple airfoils is also continuing. Finally, 
the search for Improvements in the iterative solutions involved continues in 
arder to improve the computational efficiency of all the solutions. 

Symbols are defined in an appeiidix. 


BOUNDARY-FITTED COORDINATE SYSTEM 


General Formulation 


The basic idea of the boundary-fitted coordinate systems is to numeri- 
cally generate a curvilinear coordinate system having some coordinate line 
coincident with each boundary of the physical region of interest, regardless 
of the shape of these boundaries. This is done by taking the curvilinear 
coordinates to be solutions of an elliptic partial differential system, with 
constant values of one of the curvilinear coordinates specified as Dirichlet 
boundary conditions on each boundary. Values of the other coordinate are 
either specified in a mono tonic variation over a boundary as Dirichlet bound- 
ary conditions, or are determined by Neumann boundary conditions thereon. In 
the latter case, the curvilinear coordinate lines can be made to Intersect the 
boundary according to some specified condition, such as normalcy or parallel 
to some given direction. It is also possible to exercise control over the 
spacing of the curvilinear coordinate lines in the field in order to concen- 
trate lines in regions of expected high gradients. 

In any case, the numerical generation of the coordinate system is done 
automatically for any shape boundaries, requiring only the input of points on 
the boundary. The technique has been described in detail in earlier reports 
(ref. 1 and 7) and the computer code, together with instructions for and ex- 
amples of its use in the numerical solution of partial differential equations, 
is given in reference 8. 

As mentioned previously, the curvilinear coordinates are generated by 
solving an elliptic system of suitable form. One such system is 


c xx + Sy " p(c ’ n) 


(la) 


n xx + n yy * ««•»»> 


(lb) 


with Dirichlet boundary conditions, one coordinate being specified to be equal 
to a constant on the body and equal to another constant on the outer boundary, 
with the other coordinate varying monotonically over the same range around 
both the body and the outer boundary. 
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Since it is desired to perform all numerical computations in the uniform 
rectangular transformed plane, the dependent and independent variables must be 
Interchanged in Eq. (1). This results in the coupled system 


ox. r - 2Bx r + yx 

€5 nn 


where 


- J 2 [x^P(C,n) + x^QU.n)] 


(2a) 


ay « - 26y en + Yy nn 


- t2 


J [y^P(C.n) + y n Q(C,n)J 


(2b) 


a * x 2 + y 2 
n 


T ■ + y\ 


6 - x 5 x„ + y { y„ 


J * Vn - Ve 


(Transformation relations are given in ref. 8). 


The system described by eq. (2) is a quaslllnear elliptic system for the 
coordinate functions x(C,n) and y(£,n) in the transformed plane. This set is 
considerably more complex than the linear system specified by eq. (1), but the 
boundary conditions are specified on straight boundaries, and the coordinate 
spacing in the transformed plane is uniform. 


The coordinate lines may be spaced as desired around the boundaries, 
since the assignment of the coordinate values to the [x,y] boundary points is 
arbitrary. Control of the radial spacing of the coordinate lines is accom- 
plished by varying the functions P(£,n) and Q(C,n) in equations (2). 


Automatic Concentration of Coordinate Lines into a Boundary Layer 

Consider the coordinate system generation equations (2) applied to the 
one-dimensional case of straight boundaries parallel to the x-axis. With 
n ■ constant on these boundaries, and the C-lines being normal to the bound- 
aries, we have y ^ ■ y - y ■ 0, and the x-equation is identically zero so 
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that the coordinate equations reduce to 


yy 


nn 


+ ■ 


( 3 ) 


P* 1 * ® an ba made a P«rfect differential by chooaing 
trol function Q to be 8 


the fora of the con- 


Q(n) i 



( 4 ) 


where the minus sign has been Introduced 
can be integrated to yield 


merely for convenience. 


Then eq. (3) 


y(n) - Cjf(n) + c 2 


(5) 


The constants of integration may be evaluated from 
y(l) * y v y(J) - yj Sc that 


the boundary conditions: 



This equation can then be solved for f(n) to yield 


f(n) - f(l) y(n) ~ y i 
f(J) - f(l) " yj - y x 


( 6 ) 


(7) 


JfS • 2 S arblt 5 ary d ^ltion of f(l) and f(J), will yield the required 
ti l* he requlred Q (r »> via substitution in eq. (4), to produce a 

evaluation of f(n), however, by solving eq. (3) for Q to produce 


Q(n) 



( 8 ) 


rithmS fnn^M u 8m °c tb fu 2 Ctl ° nS for y( *>’ 8uch a8 exponentials, loga- 
tr^iin ?? * ' hy P" bollc functions, etc., may be found which will concen- 
trate lines near y x with a spread out to y^ However, since the boundary 
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layer thickness at high Reynolds number is only a very small fraction of the 
distance to outer boundary of the computational field, such smooth functions 
cannot allow the lines to spread rapidly enough outside of the boundary layer. 
The result is that nearly all of the lines in the field will be within a few 
boundary layer thicknesses of the body, with a great gap near the outer 
boundary. 

Therefore, a composite function was used for y(n) , formed by Joining a 
logarithmic function to a quart ic polynomial near the edge of the boundary 
layer. This function is constructed as follows: assume that it is desired 

to space the lines in the boundary layer such that the change in velocity 
from each to the next is the same. Let the velocity profile in the boundary 
layer be approximated by the exponential 


u(y) * 1 - e~ cy 


(9) 


Let the edge of the boundary layer be defined by 


u * 0.99 at y * S 
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| 

i 

4 
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yk 
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Then the decay factor c will be given by 

c = — | In (0.01) 


Now solve eq. (9) for y(u): 

y(u) ■ - ^ In (1-u) ( 10 ) 

In order to achieve the same velocity change from each line to the next, take 

u “ 0.99 (— _^) where Is the line at the edge of the boundary layer. Sub- 
stitution in eq. (10) then yields 

y(n) » - £ In [1 - 0.99 (^-)],1 < n < n 6 (11) 

6 


Let this logarithmic function be joined to a quartic polynomial at some 
line inside or at the outer edge of the boundary layer. Thus with the function 
at n*N, the polynomial is of the form 


I 
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(12) 


y<n) “ y'(N) [n-N] + | y"(N) [n-N] 2 

+ i y"(N) [n-N] 3 + a (n-N) 4 + y(N) 


N < n < J 


Here y'(N) Is functional notation, etc. The derivatives are determined by 
differentiation of eq. (11), with evaluation at n*N. The remaining coefficient 
"a" is used to satisfy the boundary condition at the outer boundary, y(J) ■ y_. 
Thus J 

yj - y(N) - y'(N) [J-N] - j y"(N) [J-N] 2 - i y"(N) [J-N] 3 

(J-N) 4 


Note that the junction to the polynomial need not occur at the edge of the 
boundar layer, but anywhere inside it. It has been found advantageous to 
place the junction two or three lines Inside the boundary layer. 

Thus if the boundary layer thickness, 6, and the number of lines therein, 
rig, are specified, along with the distance to the outer boundary, y^, and the 

total number of lines J, and the junction line N, the control function Q(n) can 
be evaluated from 


Q(n) m - jZ 


0.99 

V 1 


1 - 0.99 


(J=V>’ 

V 1 


1,2 , N <_ r)g 


(13a) 


Q(n) •-£ ^ 

J y’(N) + y" 


y 1 (N) ± y"(N) [n-N] + 12a[n-Nl 2 


y’(N) + y"(N) [n-N] + \ y"(N) [n-N] 2 + 4a[n-N] 3 


n - N, N + 1, , J (13b) 


with the required derivatives given by 


y 


(k) 


(N) 


(k-1) I , 0.99x k 
c _ N -1^ 

[1 - 0.99 (~)] k ’ 


1,2,3 


(14) 


< 
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and y(N) by 


y(N) - - i ln[l - 0.99 (^~j-)] (15) 

6 


Although this analysis is developed for the one-dimensional case of a flat 
boundary, the same general results will be achieved by its use with curved 
boundaries since curvature tends to affect both the boundary layer thickness 
and the line control in the same way. Thus convex curvature thins the bound- 
ary layer but also causes the lines to concentrate to a greater degree near the 
boundary. 

An example of coordinate systems generated with this concentration of lines 
ir boundary layers is shown in figure 1. 

The exponential velocity distribution in eq. (9) may be replaced by the 
Blasius distribution, and this has been done in practice since a more even 
spacing is obtained in the lower part of the boundary layer in that case. 


Truncation Error Induced by the Coordinate System 

Some attention must be paid to the rapidity of the change of coordinate 
line spacing with strong attraction, else truncation error in the form of 
artificial diffusion may be introduced as follows: Consider the finite dif- 

ference approximation of a first derivative with variation only in the x- 
direction to which the £-lines are normal. Then 


f 

x 


Vi 

x.y 

5 7 n 


ft 

x « 


The difference approximation then would be 


f 

x 




+ T, 


(16) 


(17) 


where is the local truncation error. Taylor series expansions of f^ + ^ and 
^1 1 a ^ out c ^ en yield* after some algebraic rearrangement. 


T «**>! <X i+l + X l-1 - 2x i> 


(18) 
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But the last factor 


is simply the difference 


approximation of x^ 


so that 


T 


2 x ce 


f 

XX 


This truncation error thus introduces a numerical diffusive effect in the 
difference approximation of first derivatives. Care must therefore be taken 
that the second derivatives of the physical coordinates (i.e., the rate of 
change of the physical spacing between curvilinear coordinate lines) are not 
too large in regions where the dependent variables have significant second 
derivatives in the direction normal to the closely spaced coordinate lines. 

Just what is a permissible upper limit to the rate of change of the line 
spacing is problem dependent. Consider, for instance, viscous flow past an 
infinite flat plate parallel to the x-direction. Here the velocity parallel 
to the wall changes rapidly from zero at the wall to its free stream value 
over a small distance that is of the order of 1 , where R is the Reynolds 

u „ x nr 

number, R ■ - , based on freestream velocity, U w , the distance from the lead- 

ing edge of tne plate, x, and the kinematic viscosity, v. The equation for 
the time rate of change of the velocity parallel to the wall is 

u„ ■ -uu - vu + —■ (u + u ) (19) 

t x y R xx yy v ' 


Recalling that the large spacial variation in velocity occurs in the y-direction, 
coordinate lines would be contracted near the plate. The truncation error 
introduced by this contraction would be 


T(v) 


( - 2 y nn )u yy 


( 20 ) 


This introduces a negative numerical 
both positive. 


viscosity (-ly ) , since v and y 
2 nn 7 nn 


are 


The effective viscosity is thus reduced (effective Reynolds number in- 
creased), so that the velocity gradient near the wall is steepened. Therefore 
care should be taken that y is limited so that the numerical viscosity 

v nn 1 

(- '2 * 8 not significant in comparison with the physical term (— ). The 


situation is mitigated somewhat of the fact that the numerical viscosity is 
proportional to the small velocity normal to the wall, this velocity being of 
order 1 . Actually this limit is conservative, since the normal velocity 

drops to zero at the wall and only attains the order -4 on the outer portion 

/IT" 


190 


of the region of large gradient of velocity parallel to the wall, where u is 
very snail. W 


Transformation of Infinity to an Interior Circuit 

With multiple-bodies involving two closely-spaced bodies, the coordinate 
line spacing in the area between the bodies will not be too good with the type 
of coordinate system generated by the original TOMCAT code (ref. 8) as illus- 
trated in figure 2. This results from the smoothness of solutions of Poisson 
equations, the contours of which tend to avoid concave regions. Coordinate 
system control can improve the situation somewhat as shown in reference 3. 

However, with two bodies, an even better system can be obtained by first 
transforming the remote outer boundary circle to a small interior circle, or 
even a point, by the analytical complex transformation Z" = — . Here 
Z = x + iy is a point in the physical plane, while Z' *■- x' + Z iy' is a cor- 
responding point in the transformed plane. These original and transformed 
boundary circuits are illustrated in figs. 3&4 for an airfoil-flap combination. 

The coordinates (x',y') of the transformed boundary circuits are then in- 
put to the TOMCAT code in the usual manner, with the body circuits comprising 
the entire top and bottom sides of the rectangular transformed plane and one 
half of the small outer boundary circle appearing on a portion of the left and 
right sides, pie result of the TOMCAT code, in the form of E, and n lines 
drawn in the x' - y" plane, is shown in fig. 5 for a wing-slat. 

The inverse analytic transformation, Z * jt , back to the x - y plane then 
produces coordinate lines of constant E, and n in the x - y plane as shown in 
figs. 6-8, the second figure being an expanded view of the slot region. A simi- 
lar coordinate system for a leading edge slat configuration is shown in figs. 
9&10. 


INCOMPRESSIBLE SOLUTION 


liie incompressible solution is written in the primitive variable formula- 
tion, and some description has been given in reference 3. The solution is 
completely implicit, so that all equations are solved simultaneously at each 
time step. The construction essentially parallels that described for the com- 
pressible solution in the next section, with the exception that the Poisson 
equation for the pressure has no time derivative, of course. On the airfoil 
surface, the pressure is determined by iteratively adjusting the pressure at 
each point on the surface in proportion to the divergence of the velocity at 
the same point, so that upon convergence the continuity equation is s a tisfied 
on the airfoil surface (ref. 3). 

This code has been written with great generality, so that it can serve as 
a research tool with which different numerical representations can beevaluated. 
The following features can be selected simply by input options without changing 
the code: 
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(a) 2-point central convective derivatives 

(b) 6-point central convective derivatives 

(c) 2-point upwind convective derivatives 

(d) 3-point upwind convective derivatives 

(e) 4-point cross derivatives 

(f) 7-point cross derivatives 

(g) product-of-average nonlinear terms 

(h) average-of-product nonlinear terihs 

(i) fully nonlinear 

(j) second-order linearization 

(k) conservative form 

(l) nonconservative form 

(m) time-dependent form 

(n) steady-state form 

Input options can also select several starting procedures as follows: 

(a) Impulsive start from rest 

(b) Impulsive start from potential flow 

(c) Accelerating start from rest 

(d) Decreasing penetration start 

(e) Increasing Reynolds number start 

The pressure boundary condition can be selected as from the continuity 
equation as discussed above or from the normal momentum equation applied on the 
body. These boundary conditions can be either first or second order. The time 
derivative can be either first or second order, and first or second-order pro- 
jection can be used for the initial guess for the iteration at each time step 
if desired. Diffusion based on the change between time steps can be activated 
if desired, as can be flux-corrected transport or a two cell wave length filter. 

The code also contains the following choices of iteration scheme, with 
selection by input option: 

(a) point SOR with computed acceleration parameters 

(b) line SOR, either rows or columns 

(c) ADI 

(d) Stone strongly implicit 

The point SOR iteration uses a variable acceleration parameter field, the 
parameter at each point being continually adjusted to be equal to the locally 
linearized optimum. This feature is important at high Reynolds number since 
the optimum acceleration parameter depends on both the local velocity and the 
mesh spacing. Finally, a multi-grid iteration form of this same code has been 
written, and three variations of the multi-grid algorithm are now under study. 

These codes treat the ful 1 incompressible Navier-Stokes equations for any 
number of multiple bodies and thus can handle all configurations produced by 
the TOMCai coordinate code. Experimentation is now in progress with all these 
options at high Reynolds number. 
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The incompressible Navier-Stokes solution ha* also been written in the 

f0rmulatlon 6) and in the in t£r “differential 

(vorticity-velocity) formulation (ref. 9), as well as the primitive va£!ble 

"fr lts t° r s , la t about haca 8in * ie a »d 

hi!h 11 u C high Reynolds number (10 6 ) have been obtained. These 

i™? ^e 1 ? 8 n r b6r result8 » however, predict drag that is jout twice the 
xperimental value, indicating the need for further development. 

for a u eBX1 \ tS ^S descrl hed inref. 9 are shown in figures 11-13 

Mole numb * r 0f /° ‘ Som f lower Reynolds number results for a mul- 

tiple airfoil are shown in figs. 14-16. 


COMPRESSIBLE SOLUTION 


Numerical Formulation 

com P re88lb l e solution is based on tho full Navier-Stokes equations in 

5 0r V n ,.. the tran8forraed coordinates. At the interior field 
rtf H?ff second-order backward time, central space scheme is used to represent 
he differential equations. Along the body surface, the continuity equation is 
represented using second-order, one-sided space differences. Q ” is 


Turbulence Model 


The compressible code 
viscosity model as used in 


also includes an algebraic two-layer turbulent eddy 
reference 10 (termed Model 3 therein) . 


ntqu.L CO 


The compressible code is still under development, but preliminary results 
are shown in figs. 17-21 for a NACA 0018 airfoil at 5° anele of attack 
number of 0.9, Reynolds number of 20,000, and surface temperature of 0.8* 

The Mach contours in figure 18 show a rather diffuse shock and the thick 
boundary layer resulting from the somewhat low Reynolds number. The velocity 
vectors of figure .19 indicate separation on the upper surface. The peak in 7 
the density -profiles ahead of and behind the airfoil (fig, 21 ) are compression 

Jotential 8 fi n WSVe8 resulting from the impulsive start from incompressible 


CONCLUDING REMARKS 


„ ™ ? 0th Ji he lnco "*P r f 88ibl e and compressible solutions are still under devel- 
pment. The search also continues for the most appropriate multi-element 
coordinate configurations. As noted, a wide variety of numerical procedures 

n^ d r 6 floT 1Uatl0n ^ ° rder C ° d6Vel0P 80 8fflclent treatment of Sigh Reynolds 
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SYMBOLS 

Curvilinear Coordinate 
Curvilinear Coordinate 
Cartesian Coordinate 
Cartesian Coordinate 
Coordinate Control Function for £ 

Coordinate Control Function for rt 
Metric Function Defined with equations (2) 

Metric Function Lefined with Equations (2) 

Metric Function Defined with Equations (2) 

Jacobian, Also Total Number of Coordinate Lines Around Body 

General Function 

Constants of Integration 

Exponential Decay Factor 

Velocity Component 

Boundary Layer Thickness 

Free-stream Mach Number 

Junction Line 

Coefficient in Quartic Polynomial 

Truncation Error 

Time 

Reynolds Number 
Velocity Component 
Complex Variable 
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Subscripts 

K Denotes Differentiation with Respect to £ 

D Denotes Differentiation with Respect to ti 

x Denotes Differentiation with Respect to x 

y Denotes Differentiation with Respect to y 

t Denotes Differentiation with Respect to t 

6 Denotes Value at Edge of Boundary Layer 

J Denotes Value at Outer Boundary 

N Denotes Value at Junction Line 
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Figure 2.- Wing-flap coordinate 


system. Poor spacing in concave region. 
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TRAILING EDGE DETAILS 


Figure 13.- Velocity vector field of NACA 0018 airfoil. Leading- and 
t railing-edge details; R - 1,000,000; 0 - 0°; T - 2.25. 



R « 1000 
T • 3.02 











Figure 17.- Coordinate system of NACA 0018 airfoil. R - 20,000 

CONTOUR INTERVAL: 0.05 



Figure 18.- Mach contours. M w - 0.9: R - 20,000; 0 - 5°. T - 1 








